% Basic Introduction:
% This function computes the transition prob. matrix from an initial grid X
% to a given a fixed histgram approximation grid HistGrid. In this function,
% Xp is the transition destination of X. 
% Function form:
% [varargout]=ApproximateTrProbMat_UnitDim(HistGrid,Xp,TypeFlag,OutputFlag)
% Input Variables:
%       HistGridCell        Histgram grid on each dimension
%       XpMat               Transition destination from X
% Special Notes:
% This function use the external Mex function lookup from CompEcon toolbox.
function [varargout]=DistApp_Hist_TrMatApprox(HistGridCell,XpMat,OutputFlag)
    if nargin<3
        OutputFlag  =   'TrMat';
    end
    switch iscell(HistGridCell)
        % Multi-dimensional Case
        case 1
            Dim     =   length(HistGridCell);
            if Dim==size(XpMat,2)
                UnitTrMatCell   =   cell(Dim,1);
                UnitNum         =   zeros(Dim,1);
                for ii=1:Dim
                    UnitNum(ii)     =   length(HistGridCell{ii});
                    UnitTrMatCell{ii} ...
                                    =   ApproximateTrProbMat_1d(HistGridCell{ii},XpMat(:,ii));
                end
            else
                error('Dimension Mismatch!');
            end
        % 1-dimension case
        case 0
            Dim     =   1;
            if (size(HistGridCell,2)==1) && (size(XpMat,2)==1)
                UnitNum         =   length(HistGridCell);
                UnitTrMatCell   =   ApproximateTrProbMat_1d(HistGridCell,XpMat);
            else
                error('Dimension Mismatch!');
            end
    end
    switch OutputFlag
        case 'UnitTrMatCell'
            varargout{1}    =   UnitTrMatCell;
        case 'TrMat'
            if Dim==1
                varargout{1}    =   UnitTrMatCell;
            else
                varargout{1}    =   DistApp_Hist_TrMatAssemble(UnitTrMatCell,UnitNum);
            end
    end
end

function [UnitTrMat]=ApproximateTrProbMat_1d(HistGrid,XpVec)
%% Preliminary
L_Grid      =   length(HistGrid);
L_Xp        =   length(XpVec);
%% Index
Ind         =   celookup(HistGrid,XpVec);

I_leftout   =   (Ind==0);
I_rightout  =   (Ind==L_Grid);
% If the destination is to the left of left histgram boundary, 
% its destination interval left end will be assigned as the left boundary

Ind_left    =   Ind;
Ind_left(I_leftout) ...
            =   1;
% If the destination is to the right of right histgram boundary, its
% destination interval right end will be assigned as the right boundary

Ind_right  =   Ind+1;
Ind_right(I_rightout) ...
            =   L_Grid;
% For those with destination point outside of boundary, their left/right
% destination interval end points are the same.

%% Probability
Gap_left    =   XpVec-HistGrid(Ind_left);
Gap_right   =   HistGrid(Ind_right)-XpVec;
TrProb_left =   Gap_right./(Gap_left+Gap_right);
TrProb_left(I_leftout) ...
            =   0;
TrProb_left(I_rightout) ...
            =   1;
TrProb_right=   1-TrProb_left;

TrProbInd   =   [[Ind_left;Ind_right],[(1:1:L_Xp)';(1:1:L_Xp)']];
TrProb      =   [TrProb_left;TrProb_right];

%% Transition Probability Matrix
UnitTrMat   =   sparse(TrProbInd(:,1),TrProbInd(:,2),TrProb,L_Grid,L_Xp);
end